/*g_T8_col_joint_pvals.do ******************************************************

For Table 8 (Subgroup analysis), find the p-values for a test
where the null hypothesis is that there are no differences in estimates across
outcomes for a given subgroup

*******************************************************************************/

estimates clear

* Set switches
local mcas_stacked		1
local frpl_group  		1
local boygirl_group 	1
local race_group 		1

* Create dataset for stacked MCAS file
if `mcas_stacked' == 1 {
	use "$stata_data_analysis/full_samp.dta", clear

	keep studentno att in_sims in_K1 soffany yr_asgn hispanic asian female black white ///
	age_PreK_enr num_ranked first_walk pop med_inc_fam pov share_black share_white share_hispanic ///
	std_scaled_e3-std_scaled_e8 ///
	std_scaled_e10 std_scaled_m3-std_scaled_m8 std_scaled_m10 ELC soff* pform* ///
	frpl not_frpl nLC sp_biling inc1 inc2 inc3 inc4  ///
	std_e3-std_e8 std_e10 std_m3-std_m8 std_m10

	reshape long std_scaled_e std_scaled_m behav_index, i(studentno) j(grade)

	gen model = 12

	gen male = female == 0

	save "$stata_data_analysis/mcas_stacked.dta", replace
}

* Subgroup analysis by free/reduced lunch status
if `frpl_group' == 1 {
	use "$stata_data_analysis/final_sample.dta", clear
	drop grade

	gen male = female == 0

	local outcome_list any_6m_asgn any4yr_nsc_6m_asgn grad_ontime_asgn  ///
					   everany_nsc everany4yr_nsc evergrad_nsc has_sat_took_10 ///
					   grad_hs hs_ever_anysusp hs_ever_incarcerated ///
					   behav_index_tot  std_scaled_m

	local group_var "frpl not_frpl"


	matrix pval = J(60,6,.)
	local i = 1

	gen id = _n

	* Loop over outcomes
	local outcome_count = 1
	expand 11
	gen y = .
	bysort id: gen model = _n

	* Append MCAS stacked dataset:
	append using "$stata_data_analysis/mcas_stacked.dta"

	* Loop over subgroups
	foreach group of local group_var {
		foreach v of varlist att soffany pformany {
			gen `v'_`group' = `v' * `group'
		}

		egen p_`group' = group(pformany_`group')
	}

	local endog_vars ""
	local instr_vars ""

	foreach y of local outcome_list {
		gen model`outcome_count' = 1 if model == `outcome_count'
		drop if in_sims == 0 & inlist(model, 8, 9, 10, 11) // restrict to in_sims for SIMS outcomes
		display in red "outcome: `y'"
		display in red "outcome count: `outcome_count'"
		replace y = `y' if model == `outcome_count'

		foreach group of local group_var {
			display in red "group var: `group'"
			foreach v of varlist att soffany pformany {
				gen `v'_`group'_model`outcome_count' = `v'_`group' * (model == `outcome_count')
			}

			egen p_`group'_`y' = group(pformany_`group'_model`outcome_count')

			local endog_vars "`endog_vars' att_`group'_model`outcome_count'"
			local instr_vars "`instr_vars' soffany_`group'_model`outcome_count'"
			local psc_controls "`psc_controls' pformany_`group'_model`outcome_count'"
		}



		local outcome_count = `outcome_count' + 1
	}

	display in red "endog_vars: `endog_vars'"
	display in red "instr_vars: `instr_vars'"
	display in red "pscore controls: `psc_controls'"


	ivreghdfe y (`endog_vars' = `instr_vars') if in_K1 & !inlist(pformany, 0, 1), robust absorb(pformany_*_*) cluster(studentno) // dummy pscore version

	estimates store joint

	test (_b[att_frpl_model1] = _b[ att_not_frpl_model1])	///
		 (_b[att_frpl_model2] = _b[ att_not_frpl_model2])	///
		 (_b[att_frpl_model3] = _b[ att_not_frpl_model3])	///
		 (_b[att_frpl_model4] = _b[ att_not_frpl_model4])	///
		 (_b[att_frpl_model5] = _b[ att_not_frpl_model5])	///
		 (_b[att_frpl_model6] = _b[ att_not_frpl_model6])	///
		 (_b[att_frpl_model7] = _b[ att_not_frpl_model7])	///
		 (_b[att_frpl_model8] = _b[ att_not_frpl_model8])	///
		 (_b[att_frpl_model9] = _b[ att_not_frpl_model9])	///
		 (_b[att_frpl_model10] = _b[ att_not_frpl_model10])	///
		 (_b[att_frpl_model11] = _b[ att_not_frpl_model11])	///
		 (_b[att_frpl_model12] = _b[ att_not_frpl_model12])

		qui: cap estadd scalar f = r(F)
		qui: cap estadd scalar p = r(p)

		cd "$results/tables/joint_f_subgroup"
	print_results  att_frpl_model* att_not_frpl_model* , model_names( joint ) hidevars_print( 0 ) outfile("joint_frpl_notfrpl.csv") ///
			stat_names( "f" "p") stat_labs( "F test" "p-value" )
}

* Subgroup analysis by gender
if `boygirl_group' == 1 {
	use "$stata_data_analysis/final_sample.dta", clear
	drop grade

	gen male = female == 0

	local outcome_list any_6m_asgn any4yr_nsc_6m_asgn grad_ontime_asgn  ///
					   everany_nsc everany4yr_nsc evergrad_nsc has_sat_took_10 ///
					   grad_hs hs_ever_anysusp hs_ever_incarcerated ///
					   behav_index_tot  std_scaled_m

	local group_var "male female"

	matrix pval = J(60,6,.)
	local i = 1

	gen id = _n

	* Loop over outcomes
	local outcome_count = 1
	expand 11
	gen y = .
	bysort id: gen model = _n

	* Append MCAS stacked dataset:
	append using "$stata_data_analysis/mcas_stacked.dta"

	* Loop over subgroups
	foreach group of local group_var {
		foreach v of varlist att soffany pformany {
			gen `v'_`group' = `v' * `group'
		}

		egen p_`group' = group(pformany_`group')
	}

	local endog_vars ""
	local instr_vars ""

	foreach y of local outcome_list {
		gen model`outcome_count' = 1 if model == `outcome_count'
		drop if in_sims == 0 & inlist(model, 8, 9, 10, 11) // restrict to in_sims for SIMS outcomes
		display in red "outcome: `y'"
		display in red "outcome count: `outcome_count'"
		replace y = `y' if model == `outcome_count'

		foreach group of local group_var {
			display in red "group var: `group'"
			foreach v of varlist att soffany pformany {
				gen `v'_`group'_model`outcome_count' = `v'_`group' * (model == `outcome_count')
			}

			egen p_`group'_`y' = group(pformany_`group'_model`outcome_count')

			local endog_vars "`endog_vars' att_`group'_model`outcome_count'"
			local instr_vars "`instr_vars' soffany_`group'_model`outcome_count'"
		}



		local outcome_count = `outcome_count' + 1
	}

	display in red "endog_vars: `endog_vars'"
	display in red "instr_vars: `instr_vars'"


	ivreghdfe y (`endog_vars' = `instr_vars') if in_K1 & !inlist(pformany, 0, 1), robust absorb(pformany_*_*) cluster(studentno) // dummy p-score version


	estimates store joint

	test (_b[att_male_model1] = _b[ att_female_model1])	///
		 (_b[att_male_model2] = _b[ att_female_model2])	///
		 (_b[att_male_model3] = _b[ att_female_model3])	///
		 (_b[att_male_model4] = _b[ att_female_model4])	///
		 (_b[att_male_model5] = _b[ att_female_model5])	///
		 (_b[att_male_model6] = _b[ att_female_model6])	///
		 (_b[att_male_model7] = _b[ att_female_model7])	///
		 (_b[att_male_model8] = _b[ att_female_model8])	///
		 (_b[att_male_model9] = _b[ att_female_model9])	///
		 (_b[att_male_model10] = _b[ att_female_model10])	///
		 (_b[att_male_model11] = _b[ att_female_model11])	///
		 (_b[att_male_model12] = _b[ att_female_model12])

		qui: cap estadd scalar f = r(F)
		qui: cap estadd scalar p = r(p)

		cd "$results/tables/joint_f_subgroup"
	print_results  att_male_model* att_female_model* , model_names( joint ) hidevars_print( 0 ) outfile("joint_male_female.csv") ///
			stat_names( "f" "p") stat_labs( "F test" "p-value" )
}

* Subgroup analysis by race
if `race_group' == 1 {
	use "$stata_data_analysis/final_sample.dta", clear
	drop grade

	gen male = female == 0

	local outcome_list any_6m_asgn any4yr_nsc_6m_asgn grad_ontime_asgn everany_nsc everany4yr_nsc ///
					   evergrad_nsc has_sat_took_10 grad_hs hs_ever_anysusp hs_ever_incarcerated ///
					   behav_index_tot std_scaled_m

	local group_var "black white hispanic"

	matrix pval = J(60,6,.)
	local i = 1

	gen id = _n

	* Loop over outcomes
	local outcome_count = 1
	expand 11
	gen y = .
	bysort id: gen model = _n

	* Append MCAS stacked dataset:
	append using "$stata_data_analysis/mcas_stacked.dta"

	* Loop over subgroups
	foreach group of local group_var {
		foreach v of varlist att soffany pformany {
			gen `v'_`group' = `v' * `group'
		}

		egen p_`group' = group(pformany_`group')
	}

	local endog_vars ""
	local instr_vars ""

	foreach y of local outcome_list {
		gen model`outcome_count' = 1 if model == `outcome_count'
		drop if in_sims == 0 & inlist(model, 8, 9, 10, 11) // restrict to in_sims for SIMS outcomes
		display in red "outcome: `y'"
		display in red "outcome count: `outcome_count'"
		replace y = `y' if model == `outcome_count'

		foreach group of local group_var {
			display in red "group var: `group'"
			foreach v of varlist att soffany pformany {
				gen `v'_`group'_model`outcome_count' = `v'_`group' * (model == `outcome_count')
			}

			egen p_`group'_`y' = group(pformany_`group'_model`outcome_count')

			local endog_vars "`endog_vars' att_`group'_model`outcome_count'"
			local instr_vars "`instr_vars' soffany_`group'_model`outcome_count'"
		}



		local outcome_count = `outcome_count' + 1
	}

	display in red "endog_vars: `endog_vars'"
	display in red "instr_vars: `instr_vars'"


	ivreghdfe y (`endog_vars' = `instr_vars') if in_K1 & !inlist(pformany, 0, 1), robust absorb(pformany_*_*) cluster(studentno) // dummy p-score version

	estimates store joint

	test (_b[att_black_model1] = _b[ att_white_model1] = _b[ att_hispanic_model1]) ///
		 (_b[att_black_model2] = _b[ att_white_model2] = _b[ att_hispanic_model2]) ///
		 (_b[att_black_model3] = _b[ att_white_model3]	= _b[ att_hispanic_model3]) ///
		 (_b[att_black_model4] = _b[ att_white_model4]	= _b[ att_hispanic_model4]) ///
		 (_b[att_black_model5] = _b[ att_white_model5]	= _b[ att_hispanic_model5]) ///
		 (_b[att_black_model6] = _b[ att_white_model6]	= _b[ att_hispanic_model6]) ///
		 (_b[att_black_model7] = _b[ att_white_model7]	= _b[ att_hispanic_model7]) ///
		 (_b[att_black_model8] = _b[ att_white_model8]	= _b[ att_hispanic_model8]) ///
		 (_b[att_black_model9] = _b[ att_white_model9]	= _b[ att_hispanic_model9])  ///
		 (_b[att_black_model10] = _b[ att_white_model10] = _b[ att_hispanic_model10]) ///
		 (_b[att_black_model11] = _b[ att_white_model11] = _b[ att_hispanic_model11]) ///
		  (_b[att_black_model12] = _b[ att_white_model12] = _b[ att_hispanic_model12])

		qui: cap estadd scalar f = r(F)
		qui: cap estadd scalar p = r(p)

		cd "$results/tables/joint_f_subgroup"
	print_results  att_black_model* att_white_model* att_hispanic_model*, model_names( joint ) hidevars_print( 0 ) outfile("joint_black_white_hispanic.csv") ///
			stat_names( "f" "p") stat_labs( "F test" "p-value" )
}
